{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Transport equation with source term\n",
    "\n",
    "$$\n",
    "\\renewcommand{\\DdQq}[2]{{\\mathrm D}_{#1}{\\mathrm Q}_{#2}}\n",
    "\\renewcommand{\\drondt}{\\partial_t}\n",
    "\\renewcommand{\\drondx}{\\partial_x}\n",
    "\\renewcommand{\\dx}{\\Delta x}\n",
    "\\renewcommand{\\dt}{\\Delta t}\n",
    "\\renewcommand{\\grandO}{{\\mathcal O}}\n",
    "\\renewcommand{\\density}[2]{\\,f_{#1}^{#2}}\n",
    "\\renewcommand{\\fk}[1]{\\density{#1}{\\vphantom{\\star}}}\n",
    "\\renewcommand{\\fks}[1]{\\density{#1}{\\star}}\n",
    "\\renewcommand{\\moment}[2]{\\,m_{#1}^{#2}}\n",
    "\\renewcommand{\\mk}[1]{\\moment{#1}{\\vphantom{\\star}}}\n",
    "\\renewcommand{\\mke}[1]{\\moment{#1}{e}}\n",
    "\\renewcommand{\\mks}[1]{\\moment{#1}{\\star}}\n",
    "$$\n",
    "\n",
    "In this tutorial, we propose to add a source term in the advection equation. The problem reads\n",
    "$$\\drondt u + c \\drondx u = S(t, x, u), \\quad t>0, , \\quad x\\in(0, 1),$$\n",
    "\n",
    "where $c$ is a constant scalar (typically $c=1$).\n",
    "Additional boundary and initial conditions will be given in the following.\n",
    "$S$ is the source term that can depend on the time $t$, the space $x$ and the solution $u$.\n",
    "\n",
    "In order to simulate this problem, we use the $\\DdQq{1}{2}$ scheme and we add an additional `key:value` in the dictionary for the source term. We deal with two examples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A friction term\n",
    "\n",
    "In this example, we takes $S(t, x, u) = -\\alpha u$ where $\\alpha$ is a positive constant. \n",
    "The dictionary of the simulation then reads:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import sympy as sp\n",
    "import numpy as np\n",
    "import pylbm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xt8VPW9//vXZyaZJJAECYRrEFBQbgJiCOAFQZSK3YJV2+J2W+3x2mr3Pqfd/ZVf293Tg56622r9/dq6H5Zaq931AtqjsFu3WFHUIljC/RoIl0BACASIIRCSyXzOH2smDJD7rJk1k3yej8eYuaxZ65NlmPd8v9+1vktUFWOMMcbndQHGGGOSgwWCMcYYwALBGGNMmAWCMcYYwALBGGNMmAWCMcYYwKVAEJGbRaREREpFZF4zy3xFRLaKyBYRecWN7RpjjHGPxHoegoj4gR3ATUA5sBq4S1W3Ri0zHFgE3KCqx0Wkj6pWxLRhY4wxrnKjhVAElKrqblWtA14D5py3zIPAs6p6HMDCwBhjkk+aC+sYCOyPelwOTDpvmcsARGQF4Ad+rKrvnL8iEXkIeAige/fuV40YMcKF8owxputYs2bNUVXN78h73QgEaeK58/uh0oDhwDSgAPhYRMao6olz3qS6AFgAUFhYqMXFxS6UZ4wxXYeIlHX0vW50GZUDg6IeFwAHm1hmsarWq+oeoAQnIIwxxiQJNwJhNTBcRIaKSACYCyw5b5m3gOkAItIbpwtptwvbNsYY45KYA0FVg8BjwFJgG7BIVbeIyHwRmR1ebClQKSJbgQ+A76pqZazbNsYY456YDzuNl6bGEOrr6ykvL6e2ttajqryTmZlJQUEB6enpXpdijEliIrJGVQs78l43BpUTpry8nJycHIYMGYJIU2PZnZOqUllZSXl5OUOHDvW6HGNMJ5VSU1fU1tbSq1evLhUGACJCr169umTLyBiTOCkVCECXC4OIrvp7G2MSJ+UCwRhjTHxYILTT1Vdf3eoyDzzwAFu3OlM5/eQnP2n3+7OzsztWnDHGxCCljjLatm0bI0eO9KiijsnOzubkyZOuvCcVf39jTGLFcpSRtRDaKfLtffny5UybNo0777yTESNGcPfddxMJ12nTplFcXMy8efM4ffo048eP5+677z7n/SdPnmTGjBlMmDCBK664gsWLF3vzCxljTFhKHXYaLV6DrO1pMa1bt44tW7YwYMAArrnmGlasWMG1117b+Pq///u/8+tf/5r169df8N7MzEzefPNNcnNzOXr0KJMnT2b27Nk2eGyM8Yy1EGJQVFREQUEBPp+P8ePHs3fv3ja/V1X5/ve/z9ixY7nxxhs5cOAAhw8fjl+xxhjTipRtISTD2EdGRkbjfb/fTzAYbPN7X375ZY4cOcKaNWtIT09nyJAhdp6BMcZT1kKIs/T0dOrr6y94vqqqij59+pCens4HH3xAWVmHZ6w1xhhXWCDE2UMPPcTYsWMbB5Uj7r77boqLiyksLOTll1/GLgZkjPGaHXaaQrr672+MaZ0ddmqMMSZmFgjGGGMACwRjjDFhFgjGGGMACwRjjDFhFgjGGGMACwRjjDFhFgjGGGMAC4QO+c1vfkP//v0ZP34848aN48tf/jJ79uxp8T1/+tOfmDRpEuPGjaOwsJClS5cmqFpjjGkbVwJBRG4WkRIRKRWReU28fp+IHBGR9eHbA25s1ysbN25k/vz5rF+/ng0bNjBjxgxuv/32Zifce+WVV3jqqadYvHgxGzZs4NVXX+Xee++lvLw8wZUbY0zzYg4EEfEDzwKzgFHAXSIyqolFF6rq+PDt+Vi366VNmzYxZsyYxsePPPIIhw4dYv/+/RcsW1NTw7x581i0aBH9+vUDYPjw4UybNo1ly5YlrGZjjGmNGy2EIqBUVXerah3wGjDHhfW2SCQ+t7bYvHkzo0ePPue5rKwsjh8/fsGyr732GhMmTGDQoEHnPJ+RkcGpU6c6/Pt3Vapq+82YOHEjEAYC0V+Ny8PPne8OEdkoIm+IyKAmXkdEHhKRYhEpPnLkiAuluW///v3k5OSQm5vb+Fx9fT2fffYZAPfffz933nln42ubN29m3LhxF6xnw4YNjBgxgrfeeosHH3yQOXPm8O6778b/F0hhn376KcOHD2fo0KEcOnTI63KM6XTcCISmvlef35n+X8AQVR0LvAe81NSKVHWBqhaqamF+fn6LG1WNz601GzduvKB18Pvf/54bbriBcePG8bvf/e6c13Jzc6mrqzvnuZUrV1JTU8P111/Pbbfdxm9/+1tefPFFFi5c2HoBXVAoFOInP/kJ11xzDbt27aKiooJXXnnF67KM6XTcCIRyIPobfwFwMHoBVa1U1TPhh78FrnJhu544f/zg3Xff5cknn+Spp55qcvkvfvGLLFq0iEiLZ8eOHTzwwAO88MIL+Hxnd/8TTzzBo48+Gt/iU9SLL77ID37wAxoaGpgxYwaABYIxceDGJTRXA8NFZChwAJgL/GP0AiLSX1U/Cz+cDWxzYbue2LRpE8uXL2fZsmWoKiNHjuSdd97h8ssvb3L5oqIifvjDH3LjjTdy5swZGhoa+MMf/sCUKVMAp0983rx5zJo1iwkTJiTyV0kZf/nLXwB4+umn+cY3vkG/fv1Ys2YNJSUlze53Y0z7xdxCUNUg8BiwFOeDfpGqbhGR+SIyO7zYP4vIFhHZAPwzcF+s2/XKyy+/zIEDB1izZg1r167l5ZdfbvxQqqys5JFHHmHdunU8+eSTje/5+te/zoYNG/jwww8JBAJ079698bVf/epXvPfee7zxxhs899xzCf99kp2q8vHHHwNw6623kpWVxe233w7Aq6++6mVpxnQ6dsW0FNIVf//t27czcuRI+vXrx8GDBxER/vrXvzJz5kyGDRvGjh07kLYeHmZMF2BXTDOdVqR1cN111zV+8E+fPp2+fftSWlrK+V8ajDEdZ4FgktpHH30EwNSpUxufS0tLY+7cuYDThWeMcYcFgklq0S2EaJFzPexsb2PcY4Fgkta+ffsoKyujR48e5xzqCzB+/HgASkpKCAaDXpRnTKdjgWCSVqR1cO211+L3+895LTs7m8GDB1NfX09paakX5RnT6VggmKTVXHdRROSM8S1btiSsJmM6MwsEk7QiA8rNBcKoUc6kuhYIxrjDAiGJLF++nE8++cTrMpLCiRMn2LZtG5mZmRQWNn1IdaSFsHXr1kSWZkynZYGQRCwQztqxYwcAI0aMIBAINLmMdRkZ4y4LhA744x//SFFREePHj+fhhx+mrKyM4cOHc/ToUUKhENddd13jVNa33XYbV111FaNHj2bBggWN63jnnXeYMGEC48aNY8aMGezdu5fnnnuOZ555hvHjxzf2n3dVkYHiYcOGNbtM5KztkpIS6uvrE1KXMZ2ZG5PbeSNe0xW0MpXHtm3bWLhwIStWrCA9PZ1vfvObfPjhh3zve9/jkUceYdKkSYwaNYqZM2cC8MILL5CXl8fp06eZOHEid9xxB6FQiAcffJCPPvqIoUOHcuzYMfLy8njkkUfIzs7mX//1X+Pzu6WQnTt3Ai0HQuRIo7KyMkpLS7vctB7GuC11A8Ejy5YtY82aNUycOBGA06dP06dPH3784x/z+uuv89xzz7F+/frG5X/5y1/y5ptvAs7FdXbu3MmRI0eYOnUqQ4cOBSAvLy/xv0iSi7QQhg8f3uJyo0ePpqysjK1bt1ogGBOj1A0EjyblU1Xuvffec2YzBTh16hTl5eUAnDx5kpycHJYvX857773HypUr6datG9OmTaO2thZVtQnZWtGWFgI4gfD222+zZcsW7rjjjkSUZkynZWMI7TRjxgzeeOMNKioqADh27BhlZWV873vf4+6772b+/Pk8+OCDAFRVVdGzZ0+6devG9u3bWbVqFQBTpkzhww8/ZM+ePY3rAMjJyaG6utqD3yr5tKeFADawbIwbLBDaadSoUTzxxBPMnDmTsWPHctNNN7F3715Wr17dGAqBQIDf//733HzzzQSDQcaOHcu//du/MXnyZADy8/NZsGABt99+O+PGjeOrX/0q4Mz3/+abb3b5QeXjx49TWVlJt27d6NevX4vLRs5FsENPjYmdXQ8hhXSV33/16tUUFRUxduxYNmzY0OKyke659PR0ampqSE9PT1CVxiQnux6C6VTa2l0EzpFGQ4YMsTmNjHGBBYJJOm0dUI6wcQRj3JFygZCsXVzx1pV+7/a0EKKX2717d9xqMqYrSKlAyMzMpLKyskt9OIITBpWVlWRmZnpdSkK05SzlaIMHDwagrKwsbjUZ0xWk1HkIBQUFlJeXc+TIEa9LSbjMzEwKCgq8LiMh2ttlZIFgjDtSKhDS09Mbz+41ndOJEyc4evQoWVlZDBgwoE3vsUAwxh2udBmJyM0iUiIipSIyr4Xl7hQRFZEOHRJlOr/o7qK2ns0dHQhdrTvRGDfFHAgi4geeBWYBo4C7RGRUE8vlAP8MfBrrNk3n1d7xA3DmgurevTvV1dVUVVXFqzRjOj03WghFQKmq7lbVOuA1YE4Tyz0O/AyodWGbppNq7xFGACJi3UbGuMCNQBgI7I96XB5+rpGIXAkMUtU/t7QiEXlIRIpFpLgrDhyb9g8oR1ggGBM7NwKhqY7exo5cEfEBzwDfaW1FqrpAVQtVtTA/P9+F0kyq2bt3LwCXXHJJu9538cUXAxYIxsTCjUAoBwZFPS4ADkY9zgHGAMtFZC8wGVhiA8umKfv3O43NQYMGtbLkuayFYEzs3AiE1cBwERkqIgFgLrAk8qKqVqlqb1UdoqpDgFXAbFUtbnp1pqsKhUIcOHAAoN3nXFggGBO7mANBVYPAY8BSYBuwSFW3iMh8EZkd6/pN13HkyBHq6urIy8ujW7du7XqvBYIxsXPlxDRVfRt4+7znftTMstPc2KbpfDraXQQWCMa4IaXmMjKdW+QSpB0JhP79+5OWlkZFRQWnT592uzRjugQLBJM0Ii2EjszZ5Pf7G4Mksh5jTPtYIJikEUuXEVi3kTGxskAwScMCwRhvWSCYpBEZQ+joNN92cpoxsbFAMEnDWgjGeMsCwSSFhoaGDp+UFmGBYExsLBBMUjh8+DDBYJD8/PwOXyrUAsGY2FggmKQQ6/gBnB1DKC8vp6GhwZW6jOlKLBBMUoh1/AAgIyOD/Px8GhoaOHz4sFulGdNlWCCYpOBGIAAMHOhciiMyHmGMaTsLBJMULBCM8Z4FgkkKbowhgAWCMbGwQDBJwVoIxnjPAsEkBQsEY7xngWA8FwwGOXjQuepq5AO9oywQjOk4CwTjuUOHDhEKhejbty+BQCCmdVkgGNNxFgjGc251F4EFgjGxsEAwnnPrCCOAnj17kpmZSXV1NdXV1TGvz5iuxALBeC7ybT7W8QMAEWkMFmslGNM+FgjGc24NKEdE1hNpeRhj2sYCwXjOzRZC9HqshWBM+7gSCCJys4iUiEipiMxr4vVHRGSTiKwXkb+JyCg3tms6h8gH94ABA1xZnwWCMR0TcyCIiB94FpgFjALuauID/xVVvUJVxwM/A34R63ZN52EthBRVWwuvvw7HjnldiXGJGy2EIqBUVXerah3wGjAnegFV/TzqYXdAXdiu6QRUNW5jCBYI8aEa/s8998BXvgKTJsHOnV6XZVzgRiAMBPZHPS4PP3cOEXlURHbhtBD+2YXtmk6gqqqKU6dO0b17d3JyclxZpwVC/Bw7BqNHw08H/G944w3nydJSmDIFPvnE2+JMzNwIBGniuQtaAKr6rKpeCnwP+GGTKxJ5SESKRaT4yJEjLpRmkl10d5FIU39K7WeBED/f/jb03LaCbx/6LgAv3fASDV+4BSor4aabwPZ5SnMjEMqB6FNMC4CDLSz/GnBbUy+o6gJVLVTVwvz8fBdKM8nO7e4igP79+yMijddpNu545x149aUzLGQu6QT5hXyH+97/GvNGLnbC4NQpWLjQ6zJNDNwIhNXAcBEZKiIBYC6wJHoBERke9fCLgHU4GsD9I4wA0tPT6dOnD6FQiEOHDrm23q6suhoefhimsZwCymHECAr/+iQAL/whjfr/42Fnwdde87BKE6uYA0FVg8BjwFJgG7BIVbeIyHwRmR1e7DER2SIi64FvA/fGul3TObh9hFGEdRu56/HHYd8+uD8//F3vy19m6ox0xo93xhWWBG+B7GxYvRp27fK2WNNhrpyHoKpvq+plqnqpqv6/4ed+pKpLwvf/RVVHq+p4VZ2uqlvc2K5JffHoMopenwVC7FRh0SIAZY6EA2G2813v/vudhwv+MwvmhA8utG6jlGVnKhtPxaPLCCwQ3LRjB5SVwfSL1pNRUQ4DBsCECQD84z9CRgb89a9QccNc5w3WbZSyLBCMp6zLKPm9847z89FB4dbBrbeCz/noyMuDL33JaUX8Zs9MuOgi2LQJtm71qFoTCwsE4ynrMkp+kUCYVn1ud1FEpNvod/8ZQG+/w3lg3UYpyQLBeCYYDDYeBdS/f39X120znrrj9Gn48EMYSDm99q6Fbt3ghhvOWeaGG+Dii51upZJxX3GeXLzYg2pNrCwQjGcOHz5MKBSiT58+pKenu7puuyaCOz7+2AmFbxb8l/PEF74AmZnnLOPzOU8DvP35tZCW5nQbnTyZ4GpNrCwQjGfi1V0Uvc4DBw6galNndVSku+i2zPCdW29tcrlp05yf733SDcaNg1AIiovjX6BxlQWC8Uy8BpQBcnNz6d69O6dOnaKqqsr19XcVS5cCKMMqP3WemDq1yeWuv975+be/QahosvNg5cq412fcZYFgPBOvQ07BuZSmDSzHZv9+52Chy7uVEzh+2Dmk6JJLmlx24EAYNsw5o3lvv3AgrFqVwGqNGywQjGfi2WUUvV4LhI756CPn59dGrnbuFBZCCxMQRrqNlp2a4txZtSo8V7ZJFRYIxjPx7DKKXq8FQsesXev8vD7r786doqIWl48EwlsbL4HevaGiAvbujVt9xn0WCMYz8ewyAguEWK1b5/y8/PNwC2HixBaXbxxHWCGEJtk4QiqyQDCesS6j5KXqtBCEEHl7wkcLtRIIBQVw6aXw+efw2cU2jpCKLBCMZ6zLKHnt2QNVVXB1rx34qj93Ro3bcPJgpNvob6GocQSTMiwQjCeqq6upqqoiMzOTvLy8uGzDAqHjIt1FcwaEu4taGT+IiATCn8omOgPQ69Y5Z7aZlGCBYDwR+ZAuKChw7dKZ57NA6LjIgPI1GW0bP4i49lrn5wfFOeiYMRAMnl2ZSXoWCMYTkTmGIlNMxEO/fv3w+XxUVFRQX18ft+10RpHP8MurwkcYtTEQBg92DjA6ehROjprkPLl6dRwqNPFggWA8kYhASEtLo2/fvqgqn332Wdy209lEBpTTqaNn2XrnycLCNr1X5Gzv0s7u4507GzfGoUoTDxYIxhOJCASwbqOO+Owz5xSCKdmb8NWdgcsuc65z0EaRQFhZM9a5Y4GQMiwQjCcSHQg2DXbbRQaU/2FguHVw1VXten+kd+kv+8OBsGWLM5Zgkp4FgvFEogLBpsFuv8j4waRum5w7V1zRrvdHAuGjDT3QwYOhthZKS12s0MSLBYLxRCQQBg0aFNftWJdR+0UCYXjdZudOOwMhPx+GDIGaGqgeat1GqcQCwXjCxhCSV6TLKP9Qx1oIcHYcYXf3cCBs2OBCZSbeXAkEEblZREpEpFRE5jXx+rdFZKuIbBSRZSIy2I3tmtR06tQpKisrCQQC9O7dO67bskBon+pq51KYA9MrSKusgJwc5/qY7RTpNlp9xloIqSTmQBARP/AsMAsYBdwlIqPOW2wdUKiqY4E3gJ/Ful2TuqKnrPD54ttItUBon+3bnZ83F4S7i8aMaXHK6+ZEWgj/fXCcc8cCISW48a+xCChV1d2qWge8BsyJXkBVP1DVU+GHq4D49hOYpJao7iKwS2m219atzs/rLup4dxHAhAnOtZb/UjIMzcyEffvgxAmXqjTx4kYgDAT2Rz0uDz/XnPuB/3ZhuyZFJTIQcnJyyMnJoba2luPHj8d9e6kuEghXSGyBkJ0No0ZBXYOfmqFjnCc3bXKhQhNPbgRCU+3JJr+Kicg/AYXAz5t5/SERKRaR4iNHjrhQmklGiQwEsG6j9ogEwpDPwx/eY8Z0eF2RcYSyXBtHSBVuBEI5EH3sYAFw8PyFRORG4AfAbFU909SKVHWBqhaqamF+fr4LpZlklOhAiGxn//79rSxptmxxroFw0YEtzhMdbCHA2XGEtQ12pFGqcCMQVgPDRWSoiASAucCS6AVE5ErgNzhhUOHCNk0KS3QgDB7sHNRWVlaWkO2lqpoa54qXw/x78Z2uca5/0KtXh9cXaSG8e8gGllNFzIGgqkHgMWApsA1YpKpbRGS+iMwOL/ZzIBt4XUTWi8iSZlZnugCvAmHfvn0J2V6qKilxJrabOSD27iJwGhcZGfB2ebiVsWkThEIxVmniKc2Nlajq28Db5z33o6j7N7qxHdM5JDoQLg4fR28thJZFxg+uye3YGcrnCwTgyith1ape1PYaQGblQacJcsklsRVq4sbOVDYJdebMGSoqKvD7/fTt2zch27QWQttEAmEMsR1hFC3SbXSwZ7i1sXlzzOs08WOBYBLq4EHneIMBAwbg9/sTsk0bQ2ibSCBcXOVeIEQGljepBUIqsEAwCZXo7iJwDjsVEQ4ePGhXTmvB1q0Q4Ay5n5U4ZyePHBnzOiMthOWVFgipwALBJFSiZjmNFggEGDBgAKFQyM5FaEZtLezaBSOlBGlogGHDoFu3mNc7fDj06AF/O2GBkAosEExCedFCAOs2as2OHc4BQDf2da+7CJzpKwoLYSvh6c22bwdrpSUtCwSTUJGB3UQHgh1p1LLI+MHVOe4GAjjjCKfoTuVFlzhhsHOna+s27rJAMAm1e/duAIYOHZrQ7dqRRi2LBMIojZrl1CWRcYSSNOs2SnYWCCah9uzZA3gXCNZCaFokEAYdj08LAWBV9WjnjgVC0rJAMAmjqp4FgnUZtWzrVsiliu6V+yAz0xlUdsnAgc4sGMVnrIWQ7CwQTMIcOnSI2tpaevXqRW5ubkK3bV1Gzaurc7r1xxCe0G7UKHD5HJGiItiMBUKys0AwCeNV6wDODQS7UM65du6EYBCm93ZnDqOmTJwIJVxOg/ihtBROn3Z9GyZ2FggmYbwMhJycHHr27Mnp06exa22cKzJ+MCXb/fGDiKIiqCODfZmXOTPobdvm+jZM7CwQTMJEjjC6xKPJzSLjCNZtdK5IIIxsiF8gFBY6P9fYOEJSs0AwCeNlCwHsSKPmOIGgDDzm/iGnET17OmctbwzZ5TSTmQWCSRgLhOS0dSv05zMyao45n9wDBsRlO0VFsImoayOYpGOBYBLGq5PSIqzL6ELBoHNhnCuip7yWpi6THruJE2ED4aun2eU0k5IFgkmI+vp6ysvLEZHGb+qJZi2EC+3a5cwmcV1Pdy6K05KiItjLEKp9uXDoEBw+HLdtmY6xQDAJsW/fPkKhEAUFBQQCAU9qsEC4UOMcRlnrnTtxDITx48Gf5mN9aKzzhLUSko4FgkkIr7uL4OzRTaWlpXYuQljjHEZn1jl3rrwybtvKynLyprHbaP36uG3LdIwFgkmIyICyV4ecAvTu3Zu8vDyqq6s5dOiQZ3Ukk61bIZPT9Dm+3Tk7OY4tBHDGEdYz3nlgLYSkY4FgEsLrI4wiLr/8cgC2b9/uaR3JYutWZ0DZF2qAESOcr/FxVFRkLYRkZoFgEiIZuozgbCCUlJR4WkcyaGhwrldzJfHvLoq4/npnTqMGfGhJiU1hkWQsEExCJEsLYcSIEYAFAsCePc6lM6/pFv6mnoBAuPRS6F2QRQmXO5fq3LIl7ts0bedKIIjIzSJSIiKlIjKvidenishaEQmKyJ1ubNOklmQYQwBrIUSLnBtWlB5uIYwfH/dtisD06TaOkKxiDgQR8QPPArOAUcBdIjLqvMX2AfcBr8S6PZN6qqurOXr0KBkZGfTr18/TWiwQztq4EXw0cEnNRueJBAQCOIFg4wjJyY0WQhFQqqq7VbUOeA2YE72Aqu5V1Y1AyIXtmRRTWloKON1FPp+3vZSXXnopfr+fPXv2UFtb62ktXtu4ES6nhEDwNAweDHl5CdnuDTecbSHoemshJBM3/nUOBPZHPS4PP9duIvKQiBSLSLFNUdx5bAr3TYwePdrjSiAQCDB06FBUtTGouqqNGxM7oBwxeDAcH+S0EELrNjjTYZuk4EYgNDXxSYf+D6vqAlUtVNXC/Pz8GMsyySISCFfE+Rj3trJuIzh50pm24ipf4gMB4Iqb+nGYPvhrPoe9exO6bdM8NwKhHBgU9bgAOOjCek0nkWyBYEcaOQf3qMI13bwJhOnTYS0TnAerVyd026Z5bgTCamC4iAwVkQAwF1jiwnpNJ5FsgWAtBKe7CJTRdd4FwkqmABBasTKh2zbNizkQVDUIPAYsBbYBi1R1i4jMF5HZACIyUUTKgS8DvxERO/i4izh27BgHDx4kKyvL80NOIywQnEAYTBnZdcehVy8Y2KFhvw4bOBD2D3QCoea9TxK6bdO8NDdWoqpvA2+f99yPou6vxulKMl3M5vClEkePHo3f7/e4Gkf09BWqisRp/v9ktnEjXMMK58HkyXG7BkJL+tw6idBzQtb2dc4Zy3GeNsO0zs5UNnGVbN1FAH369KFHjx5UVVVRUVHhdTkJp+oEwnV87Dxx3XWe1PHFu3LZxBWkherR4jWe1GDOZYFg4ioZA0FEuvTAcnk5nDgB0/zhQLj2Wk/quOYaWJ/pdBsd+v+s2ygZWCCYuErGQICuPevpxo2QRyUjGrZCRgYUFnpSh98PDZOuBuDEOzawnAwsEEzcqGrjGEKyBUKknrVr13pcSeKdM34waZITCh4Z9jUnEPqUfmInqCUBCwQTN/v27ePzzz8nPz+fvn37el3OOa6+2vkg+uSTrtdVkQzjBxGT776Uo9KbXsEK9r6/29NajAWCiaNk7S4CmDBhAoFAgM2bN1NVVeV1OQlVXJw8gRDIEMoGOuG85XnrNvKaBYKJm2QOhMzMTK666ipUlU8//dTrchLmyBE4WFrDVaxBfT6YMsXrksiY5gTC6WVdr7UzvSqqAAAQNElEQVSWbCwQTNwkcyAATAl/GHalbqNVq2ASn5JOEBk/HnJzvS6JYfc4/x8uP/IxxcUeF9PFWSCYuFBVVq1aBcD4BM2z315dcRxh5crk6S6KyJxaRG16NlewmVef3Ot1OV2aBYKJi9LSUvbs2UNeXl7SBkKkhbBq1SoaGho8riYxVq6EaSx3Hnh0/sEFMjMJ3jTLub94MV3wXMGkYYFg4mLp0qUA3HTTTUkzZcX5BgwYwJAhQ6iurmZLF7i2bzAIez89zFQ+QtPTYcYMr0tqlH33bQD8Q8Nb/Pa3HhfThVkgmLiIBMIXvvAFjytpWVfqNtq4EWad/hN+QsjMmdCzp9clnXXLLYT8aUzlI179dSX19V4X1DVZIBjX1dXV8cEHHwAwc+ZMj6tpWVcKhJUrYS6vOQ+++lVviznfRRchN0zHT4irDv2ZRYu8LqhrskAwrluxYgU1NTWMGTOGgQmeVrm9ulIgbF92gGv5G8G0DJgzp/U3JJh86UsA3MZbzJsHNTUeF9QFWSAY16VKdxE4h8RmZ2eza9cudu3a5XU5cdXno9fxoZy87pakONz0ArNnAzBLllJZforHH/e4ni7IAsG4LhIIyd5dBJCWlsaXwt9Mn3/+eY+riZ+KCrixciEAOQ8kWXdRxMCBUFREpp7mZpby9NOwbZvXRXUtFgjGVYcPH2b9+vVkZmZyXZIc596ahx9+GIAXXniBuro6j6uJj3VvlTGFVZz2dcM/5x+8Lqd5X/4yAD/r+zTBoPLYYxAKeVxTF2KBYFz1xhtvAHD99deTlSJXwLr66qsZNWoUFRUVLF682Oty4iLtfz8NwO5Rt0L37h5X04KHHoJevRh2eAV35rzL++/Dd75jE6EmigWCcU1NTQ1PPPEEAPfff7/H1bSdiDS2EhYsWOBxNe4Lrt/M9Vv/gwZ8ZPw/3/e6nJbl5sL/+B8A/K7/D0lPU/7X/4KnnvK4ri7CAsG45plnnuHQoUNMnDiRO++80+ty2uWee+4hMzOT9957j9LSUq/LcY8q1ff/n6TRwGs9HuHSL431uqLWPfoo9OlD7o5i3vuX/wKcjPj1r62lEG8WCMYVR44c4Wc/+xkAP/3pT1PuwvU9e/bkq+Fj859++mmPq3HRkiX0XLuMY/Rk933zSYn/Ld27w/edlszUd3/Ir35eC8C3vgV33AGVlV4W17lZILTg2DH4+99h4UL4/e/h+efhxRfh/fdhzx7oItPftMnjjz9OdXU1s2bNYvr06V6X0yHf+ta38Pv9PPfcc/ziF7/wupzYbdyIPvooAD9iPjff3cvjgtrh4Ydh0CDYtInH/vuLLPzdSXJz4c034Yor4JlnoItdxiIxVDXmG3AzUAKUAvOaeD0DWBh+/VNgSGvrvOqqqzSRgkHVdetUf/Ur1a9+VXXgQFWngdr8LTtbdeZM1ccfV121yllHV3P8+HH92te+poCKiG7YsMHrkmLy0ksvKaCAPv/88wnb7rFjqm++qTp/vuoDD6jecovqgw+q/sd/qK5d24EVvvyyalaWKugKpujggfUaCrlednxt2qTar5/zj23yZN27ukKvuebcf3/33af6xz+qlpd7XWzyAIq1g5/lojF2yomIH9gB3ASUA6uBu1R1a9Qy3wTGquojIjIX+JKqtngw9LgxV+hTX7yf9PQQgUCI9PQQF8yR1pbao5aJ/l1Pn/bz2Wfd2Lu3O2V7u7N/fza1Z87dQEaggby8M+TlnSErMwhAKCQcPx7g2LEMPq8OnLN896w6hg2r5rLLqhg+vJrs7GCL9TRdrp5XM5w4EaCyMoujRzOoqUmj9oyfulqncScCfn+IjIwGAoEQgUADGRkhMgJnH0eeS08P0RASGoJCQ4PQEISGBp9zv/HmPFYFn0/xiSKi+HyK+MDvAwhRXV3F8eOVbNy4jqqqKtLTAtxyyxeZOLGose7oG6qIgPjAJ4A4P0WcwwpVQUN64fsAn+/sTVB8vqj36IX3BcXvdy7i7vM5PyPbCQYjP5VQAwQbIBRu6aWlQVo6rF37Ke+8818I9Qzo35fRoy+jX7+BBIPpBINpUT/91NcJdXV+RJS0NCUQCJGTW0+P3Dp69KgnK7OhyW6a+nqhvLwbe/ZkU7a3G4cOZRHS5vtzhl3yOTfMOER+7zPNLiMNDeSUlnLRpk3k7twJwIdD7uDmvf/JrC8d5dFHd7T4t5eMsg4c4MrvfpesigpCfj/Hx41nTd/pfLD1MjaWDeAMZ68H3S0rSO/eZ7joojqyshrIzGxo/JmR0eD8PfucvwWf7+zftc8HCET2/tn/X3re4+jXz33twv/HTb/eXJdds9uUc360sl3Fn5HOrJ9+c42qFja9pVZ0NEkiN2AKsDTq8f8E/ud5yywFpoTvpwFHwQmj5m6X9hnc+ld0u9nNbk3eToF+A1HYEX5qqoY/TVLuVgD6Nmh9EuzXZL+dIFeJoYWQ1qEUOddAYH/U43JgUnPLqGpQRKqAXjjB0EhEHgIeAujfoz/v+65F1Y+qj+aGO5TWR8maWkYI4ffXkpZ2krT0k6Sl1eDztXxSkrYyItfQkEldXQ/q6nsQrM9Bm6w5BBLCJyHnr119oP4Llo2u2Sd1+P21+NNq8fnqEQnikwYUCX9f8KEhH4ofDTn7S/ETCvnD+y98wxdesyISAtHwzgghOI8FPfu8Onsq/A5AGp8Tnx+/z4/P7yeQHkBFnIqjv9FEfety3uq8XzW8mvBPkbPLR77yRO/qyF975L0hxHmPnLsdCZeIQkidVk5IQUPOdiItksi3RPHJ2RYLTsshFIKG8M9Qg3KmLkRdfRA0hEhD4w0JRt1vABUUH6p+Qg2ZNDRk0hDKRNV/wf/PyM5NSztJIP04gYzjpKdXOetq6i8mFOBkzVBOnR4ICIH0Y/S4aDM+34Ut0ANZWWy46CK29+jB7sO3wvbhZGQcZsqUAD7fDU2uPxU8BSyor2fK0aOMqqqiezBI92AQv2rjMhpKoyGUSSiUfu7ffeNniIA6fyTOu8L/TzTyOOq5Nom8v43LNePCt7dWw3mvR62g1pcBwQ9aK6hZbgRCU9Wf/zu2ZRlUdQGwAKCwsFBvKP648bVQCGprnQmvamrg1KnW79fWQno6BALOgQt9+zq3YcPg4oudD4Z4OX0a1q51ZpgsLobSUudWVcXZ7z5RevRwxtCGDoWRI2HECOfn5Zcn1yzFpu1U4fhx2L/f+Xs8c8YJosGDnf/Xae3817dyJdx2mzMNxYhe8Oc/w6WXNr3s8eNw2WXO/Rdf7MvcuX+N7ZcxqSOGQ8ncCIRyYFDU4wLgYDPLlItIGtADONaejfh80K2bc8vPj6XcxMjKgmuucW4Rqs6HwuefO7e0NGe57t0hO9u7Wk18iEBennNzw5Qp8OmncOutsHkzTJrkHHXT1AwhP/4xHD0KU6cm30zXJnm58R15NTBcRIaKSACYCyw5b5klwL3h+3cC72uso9kpSAQyM6FPH6eVMmSI02KxMDBtNWQIrFgBs2Y5x+PPmOEcDh0938/f/w7PPut8ifrlL2P6wmi6mJgDQVWDwGM4A8fbgEWqukVE5ovI7PBivwN6iUgp8G1gXqzbNaarys2FJUvgX/4F6uvhwQed7qEnn4Tbb4fJk51zZL7xDRg3zutqTSqJ+bDTeCksLNTi4mKvyzAmqb34otM9VFZ29rlAAL7+dWf+H2t9dj0i0uHDTt0YQzDGeOS+++Cee+Avf3HOqB8yxJnioV8/ryszqcgCwZgU5/c7FxubPbv1ZY1pic1lZIwxBrBAMMYYE2aBYIwxBrBAMMYYE2aBYIwxBrBAMMYYE2aBYIwxBrBAMMYYE2aBYIwxBrBAMMYYE2aBYIwxBrBAMMYYE2aBYIwxBrBAMMYYE2aBYIwxBrBAMMYYE2aBYIwxBrBAMMYYE2aBYIwxBrBAMMYYExZTIIhInoj8VUR2hn/2bGa5d0TkhIj8OZbtGWOMiZ9YWwjzgGWqOhxYFn7clJ8D98S4LWOMMXEUayDMAV4K338JuK2phVR1GVAd47aMMcbEUVqM7++rqp8BqOpnItInlpWJyEPAQ+GHZ0Rkc4z1dRa9gaNeF5EkbF+cZfviLNsXZ13e0Te2Gggi8h7Qr4mXftDRjTZHVRcAC8LbLVbVQre3kYpsX5xl++Is2xdn2b44S0SKO/reVgNBVW9sYcOHRaR/uHXQH6joaCHGGGO8FesYwhLg3vD9e4HFMa7PGGOMR2INhH8HbhKRncBN4ceISKGIPB9ZSEQ+Bl4HZohIuYh8oQ3rXhBjbZ2J7YuzbF+cZfviLNsXZ3V4X4iqulmIMcaYFGVnKhtjjAEsEIwxxoR5HggicrOIlIhIqYhccKaziGSIyMLw65+KyJDEV5kYbdgX3xaRrSKyUUSWichgL+pMhNb2RdRyd4qIikinPeSwLftCRL4S/tvYIiKvJLrGRGnDv5GLReQDEVkX/ndyixd1xpuIvCAiFc2dqyWOX4b300YRmdCmFauqZzfAD+wCLgECwAZg1HnLfBN4Lnx/LrDQy5o93hfTgW7h+9/oyvsivFwO8BGwCij0um4P/y6GA+uAnuHHfbyu28N9sQD4Rvj+KGCv13XHaV9MBSYAm5t5/RbgvwEBJgOftmW9XrcQioBSVd2tqnXAazjTYUSLnh7jDZwjlSSBNSZKq/tCVT9Q1VPhh6uAggTXmCht+bsAeBz4GVCbyOISrC374kHgWVU9DqCqnfV8oLbsCwVyw/d7AAcTWF/CqOpHwLEWFpkD/EEdq4CLwueKtcjrQBgI7I96XB5+rsllVDUIVAG9ElJdYrVlX0S7H+cbQGfU6r4QkSuBQara2WfQbcvfxWXAZSKyQkRWicjNCasusdqyL34M/JOIlANvA99KTGlJp72fJ0DscxnFqqlv+ucfB9uWZTqDNv+eIvJPQCFwfVwr8k6L+0JEfMAzwH2JKshDbfm7SMPpNpqG02r8WETGqOqJONeWaG3ZF3cBL6rq0yIyBfjP8L4Ixb+8pNKhz02vWwjlwKCoxwVc2MRrXEZE0nCagS01lVJVW/YFInIjzjxSs1X1TIJqS7TW9kUOMAZYLiJ7cfpIl3TSgeW2/htZrKr1qroHKMEJiM6mLfvifmARgKquBDJxJr7ratr0eXI+rwNhNTBcRIaKSABn0HjJectET49xJ/C+hkdNOplW90W4m+Q3OGHQWfuJoZV9oapVqtpbVYeo6hCc8ZTZqtrhSb2SWFv+jbyFc8ABItIbpwtpd0KrTIy27It9wAwAERmJEwhHElplclgCfC18tNFkoErDM1O3xNMuI1UNishjwFKcIwheUNUtIjIfKFbVJcDvcJp9pTgtg7neVRw/bdwXPweygdfD4+r7VHW2Z0XHSRv3RZfQxn2xFJgpIluBBuC7qlrpXdXx0cZ98R3gtyLyf+F0kdzXGb9AisirOF2EvcPjJf83kA6gqs/hjJ/cApQCp4Cvt2m9nXBfGWOM6QCvu4yMMcYkCQsEY4wxgAWCMcaYMAsEY4wxgAWCMcaYMAsEY4wxgAWCMcaYsP8fC8C3eauzgb8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "C, ALPHA, X, u, LA = sp.symbols('C, ALPHA, X, u, LA')\n",
    "c = 0.3\n",
    "alpha = 0.5\n",
    "\n",
    "def init(x):\n",
    "    middle, width, height = 0.4, 0.1, 0.5   \n",
    "    return height/width**10 * (x%1-middle-width)**5 * \\\n",
    "                              (middle-x%1-width)**5 * (abs(x%1-middle)<=width)\n",
    "    \n",
    "def solution(t, x):\n",
    "    return init(x - c*t)*np.exp(-alpha*t)\n",
    "\n",
    "dico = {\n",
    "    'box': {'x': [0., 1.], 'label': -1},\n",
    "    'space_step': 1./128,\n",
    "    'scheme_velocity': LA,\n",
    "    'schemes': [\n",
    "        {\n",
    "            'velocities': [1,2],\n",
    "            'conserved_moments': u,\n",
    "            'polynomials': [1, LA*X],\n",
    "            'relaxation_parameters': [0., 2.],\n",
    "            'equilibrium': [u, C*u],\n",
    "            'source_terms': {u: -ALPHA*u},\n",
    "        },\n",
    "    ],\n",
    "    'init': {u: init},\n",
    "    'parameters': {\n",
    "        LA: 1., \n",
    "        C: c, \n",
    "        ALPHA: alpha\n",
    "    },\n",
    "    'generator': 'numpy',\n",
    "}\n",
    "\n",
    "sol = pylbm.Simulation(dico) # build the simulation\n",
    "viewer = pylbm.viewer.matplotlib_viewer\n",
    "fig = viewer.Fig()\n",
    "ax = fig[0]\n",
    "ax.axis(0., 1., -.1, .6)\n",
    "x = sol.domain.x\n",
    "ax.plot(x, sol.m[u], width=2, color='k', label='initial')\n",
    "while sol.t < 1:\n",
    "    sol.one_time_step()\n",
    "ax.plot(x, sol.m[u], width=2, color='b', label=r'$D_1Q_2$')\n",
    "ax.plot(x, solution(sol.t, x), width=2, color='r', label='exact')\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A source term depending on time and space\n",
    "\n",
    "If the source term $S$ depends explicitely on the time or on the space, we have to specify the corresponding variables in the dictionary through the key *parameters*. The time variable is prescribed by the key *'time'*. Moreover, sympy functions can be used to define the source term like in the following example. This example is just for testing the feature... no physical meaning in mind !"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xt8VPWd//HXJ9cJV7mEaxCCooAICBHrrVovVesWe6EtVmvtaq1a6z5+vezya/vrttaHWqvb3Z/aR0ut1bZape5W/e1isVJtrfUCKCgBwQAKAcIlYkQSIMl8fn+cGZhgAknmzJyZ5P18POZx5nLmnE8Ow7zn+/2ei7k7IiIiBVEXICIiuUGBICIigAJBREQSFAgiIgIoEEREJEGBICIiQEiBYGYXmtkaM6sxs3kdzPNZM1tlZtVm9lAY6xURkfBYuschmFkhsBY4H6gFlgCXuvuqlHkmAAuAc9x9l5kNc/ftaa1YRERCFUYLYRZQ4+7r3X0/8DBwySHzfBm4x913ASgMRERyT1EIyxgNbEp5XAuccsg8xwGY2fNAIfB9d//joQsys2uAawD69u07c+LEiSGUJyLSeyxbtmynu5d3571hBIK189yh/VBFwATgbKACeM7Mprj7u23e5D4fmA9QVVXlS5cuDaE8EZHew8ze7u57w+gyqgXGpDyuALa0M8/j7t7s7huANQQBISIiOSKMQFgCTDCzSjMrAeYCTxwyz2PARwDMbChBF9L6ENYtIiIhSTsQ3L0FuAFYBKwGFrh7tZndZGazE7MtAurNbBXwDPAtd69Pd90iIhKetHc7zZT2xhCam5upra1l7969EVUVnVgsRkVFBcXFxVGXIiI5zMyWuXtVd94bxqBy1tTW1tK/f3/GjRuHWXtj2T2Tu1NfX09tbS2VlZVRlyMiPVRenbpi7969DBkypFeFAYCZMWTIkF7ZMhKR7MmrQAB6XRgk9da/W0SyJ+8CQUREMkOB0EWnnXbaEee5+uqrWbUqOJXTLbfc0uX39+vXr3vFiYikIa/2Mlq9ejWTJk2KqKLu6devH++//34o78nHv19EsiudvYzUQuii5K/3Z599lrPPPps5c+YwceJELrvsMpLhevbZZ7N06VLmzZtHU1MT06dP57LLLmvz/vfff59zzz2XGTNmcOKJJ/L4449H8weJiCTk1W6nqTI1yNqVFtOrr75KdXU1o0aN4vTTT+f555/njDPOOPD6bbfdxt13383y5cs/8N5YLMYf/vAHBgwYwM6dO/nQhz7E7NmzNXgsIpFRCyENs2bNoqKigoKCAqZPn85bb73V6fe6O9/+9reZOnUq5513Hps3b2bbtm2ZK1ZE5AjytoWQC2MfpaWlB+4XFhbS0tLS6fc++OCD7Nixg2XLllFcXMy4ceN0nIGIREothAwrLi6mubn5A883NDQwbNgwiouLeeaZZ3j77W6fsVZEJBQKhAy75pprmDp16oFB5aTLLruMpUuXUlVVxYMPPoguBiQiUdNup3mkt//9InJk2u1URETSpkAQERFAgSAiIgkKBBERARQIIiKSoEAQERFAgSAiIgkKBBERARQI3fLzn/+ckSNHMn36dKZNm8ZnPvMZNmzYcNj3/Od//iennHIK06ZNo6qqikWLFmWpWhGRzgklEMzsQjNbY2Y1ZjavndevNLMdZrY8cbs6jPVG5bXXXuOmm25i+fLlrFixgnPPPZdPfepTHZ5w76GHHuKOO+7g8ccfZ8WKFfzud7/ji1/8IrW1tVmuXESkY2kHgpkVAvcAFwGTgUvNbHI7sz7i7tMTt3vTXW+UXn/9daZMmXLg8bXXXktdXR2bNm36wLx79uxh3rx5LFiwgBEjRgAwYcIEzj77bBYvXpy1mkVEjiSMFsIsoMbd17v7fuBh4JIQlntYZpm5dcbKlSs54YQT2jxXVlbGrl27PjDvww8/zIwZMxgzZkyb50tLS2lsbOz2399bubu2m0iGhBEIo4HUn8a1iecO9Wkze83MHjWzMe28jpldY2ZLzWzpjh07QigtfJs2baJ///4MGDDgwHPNzc1s3boVgKuuuoo5c+YceG3lypVMmzbtA8tZsWIFEydO5LHHHuPLX/4yl1xyCU899VTm/4A89tJLLzFhwgQqKyupq6uLuhyRHieMQGjvd/Whnen/Dxjn7lOBp4EH2luQu8939yp3ryovLz/sSt0zczuS11577QOtg1/96lecc845TJs2jV/+8pdtXhswYAD79+9v89wLL7zAnj17OOuss/jEJz7BL37xC+6//34eeeSRIxfQC8XjcW655RZOP/101q1bx/bt23nooYeiLkukxwkjEGqB1F/8FcCW1Bncvd7d9yUe/gKYGcJ6I3Ho+MFTTz3Frbfeyh133NHu/BdffDELFiwg2eJZu3YtV199Nffddx8FBQc3/80338xXv/rVzBafp+6//36+853v0NrayrnnngugQBDJgDAuobkEmGBmlcBmYC7w+dQZzGyku29NPJwNrA5hvZF4/fXXefbZZ1m8eDHuzqRJk/jjH//I8ccf3+78s2bN4rvf/S7nnXce+/bto7W1lV//+teceuqpQNAnPm/ePC666CJmzJiRzT8lb/zP//wPAHfeeSfXXXcdI0aMYNmyZaxZs6bD7S4iXZd2C8HdW4AbgEUEX/QL3L3azG4ys9mJ2W40s2ozWwHcCFyZ7nqj8uCDD7J582aWLVvGK6+8woMPPnjgS6m+vp5rr72WV199lVtvvfXAe770pS+xYsUK/vKXv1BSUkLfvn0PvHbXXXfx9NNP8+ijj/Kzn/0s639PrnN3nnvuOQA+/vGPU1ZWxqc+9SkAfve730VZmkiPoyum5ZHe+Pe/8cYbTJo0iREjRrBlyxbMjD/96U989KMf5dhjj2Xt2rVYZ3cPE+kFdMU06bGSrYMzzzzzwBf/Rz7yEYYPH05NTQ2H/mgQke5TIEhO++tf/wrAhz/84QPPFRUVMXfuXCDowhORcCgQJKelthBSJY/10NHeIuHJu0DI1TGPTOuNf/fGjRt5++23GThwYJtdfQGmT58OwJo1a2hpaYmiPJEeJ68CIRaLUV9f3+u+HN2d+vp6YrFY1KVkVbJ1cMYZZ1BYWNjmtX79+jF27Fiam5upqamJojyRHieM4xCypqKigtraWnL1tBaZFIvFqKioiLqMrOqouyjphBNO4O2336a6upqJEydmszSRHimvAqG4uJjKysqoy5AsSQ4odxQIkydPZuHChVRXV/PpT386m6WJ9Eh51WUkvce7777L6tWricViVFW1v0t18pxSq1atymZpIj2WAkFy0tq1awGYOHEiJSUl7c6TDITq6uqs1SXSkykQJCclB4qPPfbYDudJHrW9Zs0ampubs1KXSE+mQJCc9OabbwKHDwTtaSQSLgWC5KTkF/yECRMOO5/GEUTCo0CQnNSZFgJoHEEkTAoEyUldbSEoEETSp0CQnLNr1y7q6+vp06cPI0aMOOy8kydPBtRlJBIGBYLknNQ9jI50rQPtaSQSHgWC5JzOdhdBsKfRuHHjtKeRSAgUCJJzOjugnKRxBJFwKBAk53SlhZA63/r16zNWk0hvoECQnNOZo5RTjR07FoC33347YzWJ9AYKBMk5Xe0yUiCIhEOBIDnl3XffZefOnZSVlTFq1KhOvUeBIBKOUALBzC40szVmVmNm8w4z3xwzczNr/3zG0ut1ZZfTpNRA6G1X0xMJU9qBYGaFwD3ARcBk4FIzm9zOfP2BG4GX0l2n9FxdHT8AGDx4MH379mX37t00NDRkqjSRHi+MFsIsoMbd17v7fuBh4JJ25vshcDuwN4R1Sg/V1T2MAMxM3UYiIQgjEEYDm1Ie1yaeO8DMTgLGuPt/H25BZnaNmS01s6W98brJ0vUB5SQFgkj6wgiE9jp6D3TkmlkB8BPgG0dakLvPd/cqd68qLy8PoTTJN2+99RYA48eP79L7jj76aECBIJKOMAKhFhiT8rgC2JLyuD8wBXjWzN4CPgQ8oYFlac+mTUFjc8yYMUeYsy21EETSF0YgLAEmmFmlmZUAc4Enki+6e4O7D3X3ce4+DngRmO3uS0NYt/Qg8XiczZs3A1BRUdGl9yoQRNKXdiC4ewtwA7AIWA0scPdqM7vJzGanu3zpPXbs2MH+/fsZPHgwffr06dJ7FQgi6SsKYyHuvhBYeMhz3+tg3rPDWKf0PN3tLgIFgkgYdKSy5Iza2lqge4EwcuRIioqK2L59O01NTWGXJtIrKBAkZyRbCF0dPwAoLCw8ECTJ5YhI1ygQJGek02UE6jYSSZcCQXKGAkEkWgoEyRnJMYTudBmBDk4TSZcCQXKGWggi0VIgSE5obW3t9kFpSQoEkfQoECQnbNu2jZaWFsrLy4nFYt1ahgIhO3TJiZ5LgSA5Id3xAzg4hlBbW0tra2sodUlb770H06ZBVRWsWhV1NRI2BYLkhHTHDwBKS0spLy+ntbWVbdu2hVWapLj9dnj9dVi2DE4+GX7726grkjApECQnhBEIAKNHB5fiSI5HSHhqa+Hf/i24f8EF0NgIX/gC3HdftHVJeBQIkhMUCLnv//wfaGqCOXPgySfhttuC53/2s2jrkvAoECQnhDGGAAqETFmxAh54AIqL4dZbwQy+9jXo0weWLAGN4/cMCgTJCWoh5LYf/SjYu+j66yF5ddM+feDii4P7//Vf0dUm4VEgSE5QIOQud/jzn4P7113X9rU5c4Lpo49mtybJDAWCRK6lpYUtW4Krria/0LtLgRC+mhrYtg2GDYPjjmv72sc+BrEY/P3voE2e/xQIErm6ujri8TjDhw+npKQkrWUpEML3t78F0zPOCMYOUvXrBxddFNxXt1H+UyBI5MLqLgIFQiY891wwPeOM9l9Xt1HPoUCQyIW1hxHAoEGDiMVi7N69m927d6e9PDkYCGee2f7r//APUFISzLd9e/bqkvApECRyyV/z6Y4fAJjZgWBRKyF9dXXBGELfvjB9evvzDBgQtB7c4YUXslufhEuBIJELa0A5KbmcZMtDui85fnDqqVBU1PF8p5wSTF9+OfM1SeYoECRyYbYQUpejFkL6jjR+kDRrVjBVIOS3UALBzC40szVmVmNm89p5/Voze93MlpvZ38xschjrlZ4h+cU9atSoUJanQAhPsoXQ0fhBUjIQliyBeDyzNUnmpB0IZlYI3ANcBEwGLm3nC/8hdz/R3acDtwP/lu56pedQCyE3vfceLF8edBUlu4Q6MmoUjB4NDQ3w5pvZqU/CF0YLYRZQ4+7r3X0/8DBwSeoM7v5eysO+gC6xIQC4e8bGEBQI6XnppeDX/owZwaDykSRbCS+9lNm6JHPCCITRwKaUx7WJ59ows6+a2TqCFsKNIaxXeoCGhgYaGxvp27cv/fv3D2WZCoRwrFwZTGfM6Nz8GkfIf2EEgrXz3AdaAO5+j7sfA/wL8N12F2R2jZktNbOlO3bsCKE0yXWp3UV26GGw3aRACEfyimgnnNC5+RUI+S+MQKgFUg8xrQC2HGb+h4FPtPeCu8939yp3ryovLw+hNMl1YXcXAYwcORIzO3CdZumeZCBM7uQuIDNnBqe2WL4c9u3LXF2SOWEEwhJggplVmlkJMBd4InUGM5uQ8vBiQMNOAoS/hxFAcXExw4YNIx6PU1dXF9pyexP3rgfCwIEwcSI0NwfXT5D8k3YguHsLcAOwCFgNLHD3ajO7ycxmJ2a7wcyqzWw58HXgi+muV3qGsPcwSlK3UXrq6uDdd2HQIBg+vPPvU7dRfjvMsYed5+4LgYWHPPe9lPv/FMZ6pOfJRJdRcnmvvPKKAqGbUlsHXRnamTUruLKaAiE/6UhliVQmuoxALYR0dbW7KOnkk4PpsmXh1iPZoUCQSKnLKDd1NxBOOCFoUaxZo4HlfKRAkEhlsssIFAjdVV0dTLsaCH36wIQJ0NoKq1eHX5dklgJBItPS0nJgL6CRI0eGumyd8bT73LsfCAAnnhhMX389vJokOxQIEplt27YRj8cZNmwYxcXFoS5b10Tovh074J13guscdKfhpkDIXwoEiUymuotSl7l582bcdeqsrujuHkZJU6cG09deC68myQ4FgkQmUwPKAAMGDKBv3740NjbS0NAQ+vJ7su4OKCephZC/FAgSmUztcgrBpTQ1sNw96QbC+PHB4PKWLVBfH15dknkKBIlMJruMUperQOiadAOhoACmTAnuq5WQXxQIEplMdhmlLleB0DXJ3UUnTer+MtRtlJ8UCBKZTHYZgQKhO3bvDs5jVFoKRx/d/eUoEPKTAkEioy6j3LN+fTAdPz7o+uku7WmUnxQIEhl1GeWemppgeuyx6S0n2UJYuTK4DKfkBwWCRGL37t00NDQQi8UYPHhwRtahQOi6ZCAcc0x6yxk6FEaMgD174K230i5LskSBIJFIfklXVFSEdunMQykQum7dumCabgsB1G2UjxQIEonkOYaSp5jIhBEjRlBQUMD27dtpbm7O2Hp6krC6jEADy/lIgSCRyEYgFBUVMXz4cNydrVu3Zmw9PUmyhZBulxEoEPKRAkEikY1AAHUbdcXevbBpExQWwtix6S9PXUb5R4Egkch2IOg02Ee2YUNw6utx4yCMk89OmhSEy5tvQlNT+suTzFMgSCSyFQg6DXbnhdldBBCLBRfLicd1sZx8oUCQSCQDYcyYMRldj7qMOi/MAeUkdRvlFwWCREJjCLknrGMQUmlgOb+EEghmdqGZrTGzGjOb187rXzezVWb2mpktNrMQhqwkXzU2NlJfX09JSQlDhw7N6LoUCJ0X5jEISQqE/JJ2IJhZIXAPcBEwGbjUzA49ce6rQJW7TwUeBW5Pd72Sv1JPWVGQzglzOkGB0HnqMpIw/jfOAmrcfb277wceBi5JncHdn3H3xsTDF4HM9hNITstWdxHoUpqd1dJy8BQTlZXhLXfsWOjXD7ZtC67VLLktjEAYDWxKeVybeK4jVwFPhrBeyVPZDIT+/fvTv39/9u7dy65duzK+vny1cWMQChUVUFYW3nJ1sZz8EkYgtHcimnZ/ipnZ5UAV8OMOXr/GzJaa2dId+jnRY2UzEEDdRp2Rie6iJHUb5Y8wAqEWSN13sALYcuhMZnYe8B1gtrvva29B7j7f3avcvaq8vDyE0iQXZTsQkuvZtGnTEebsvcI+BiGVBpbzRxiBsASYYGaVZlYCzAWeSJ3BzE4Cfk4QBttDWKfksWwHwtjEeRjefvvtrKwvH2WyhaBAyB9pB4K7twA3AIuA1cACd682s5vMbHZith8D/YDfm9lyM3uig8VJLxBVIGzcuDEr68tHmTgGISn1YjmtreEvX8JTFMZC3H0hsPCQ576Xcv+8MNYjPUO2A+HoxMWB1ULoWCaOQUgaPBhGj4bNm4NLdE6YEP46JBw6Ulmyat++fWzfvp3CwkKGDx+elXWqhXB48XhmxxBA3Ub5QoEgWbVlS7C/wahRoygsLMzKOjWGcHhbtgSnvi4vhwEDMrMOBUJ+UCBIVmW7uwiC3U7NjC1btujKae3IZHdRknY9zQ8KBMmqbJ3lNFVJSQmjRo0iHo/rWIR2ZHIPoyS1EPKDAkGyKooWAqjb6HAyuYdR0sSJwcVyamqgsfHI80s0FAiSVcmB3WwHgvY06lg2uoxKS4NQcIfq6sytR9KjQJCsWr9+PQCVYZ5BrRO0p1HHstFCAHUb5QMFgmTVhg0bgOgCQS2EttyzM4YACoR8oECQrHH3yAJBXUbt27kTdu+GgQNhyJDMrkt7GuU+BYJkTV1dHXv37mXIkCEMyNQO7x1Ql1H7UruLrL3zFodILYTcp0CQrImqdQBtA0EXyjkoGwPKSUcfHRz4tmNHcMEcyT0KBMmaKAOhf//+DBo0iKamJnStjYOyNX4AQQsk2UpQt1FuUiBI1iT3MBo/fnwk60+OI6jb6KBs7WGUpG6j3KZAkKyJsoUA2tOoPdnsMoKDA8srVmRnfdI1CgTJGgVC7sl2C2HmzGD68svZWZ90jQJBsiaqg9KS1GXUVkNDsNtpWRmMHJmddU6bBiUl8MYb8O672VmndJ4CQbKiubmZ2tpazOzAL/VsUwuhrTffDKbHHAMFWfomKC2Fk04K7i9Zkp11SucpECQrNm7cSDwep6KigpKSkkhqUCC0tWZNMD3++Oyu95RTgulLL2V3vXJkCgTJiqi7i+Dg3k01NTU6FgEFgnyQAkGyIjmgHNUupwBDhw5l8ODB7N69m7q6usjqyBXJQDjuuOyuNzUQlMu5RYEgWRH1HkZJxyd+Dr/xxhuR1pEL1q4NptluIYwfD0OHBkcsv/VWdtcth6dAkKzIhS4jOBgIa5I/j3upeDy6QDCDWbOC++o2yi0KBMmKXGkhTJw4EVAgbN4cXLmsvBwGDcr++jWOkJtCCQQzu9DM1phZjZnNa+f1D5vZK2bWYmZzwlin5JdcGEMAtRCSohpQTlIg5Ka0A8HMCoF7gIuAycClZjb5kNk2AlcCD6W7Psk/u3fvZufOnZSWljJixIhIa1EgBKIOhGSX0SuvwP790dQgHxRGC2EWUOPu6919P/AwcEnqDO7+lru/BsRDWJ/kmZrE+REqKyspyNYRUB045phjKCwsZMOGDezduzfSWqIUdSAMGhTs3bRvHyxfHk0N8kFh/O8cDWxKeVybeK7LzOwaM1tqZkt1iuKe4/XEqS1POOGEiCuBkpISKisrcfcDQdUbRR0IAB/+cDBdvDi6GqStMAKhvessdWvvYnef7+5V7l5VXl6eZlmSK5KBcGLy3McRU7dRbgTC+ecH0z/9KboapK0wAqEWGJPyuALYEsJypYfItUDo7XsaNTXBxo1QWAhR7vR17rnBLqjPPx/s8STRCyMQlgATzKzSzEqAucATISxXeohcC4Te3kKoqQmOEB4/PjjzaFSGDIEZM4JB5b/+Nbo65KC0A8HdW4AbgEXAamCBu1eb2U1mNhvAzE42s1rgM8DPzaw63fVKfnjnnXfYsmULZWVlke9ymtTbAyEXuouS1G2UW4rCWIi7LwQWHvLc91LuLyHoSpJeZuXKlUAwoFxYWBhxNYHU01e4O2btDYP1XLkWCLfdpkDIFTpSWTIq17qLAIYNG8bAgQNpaGhg+/btUZeTdbkUCKefHlyg5/XXQecbjJ4CQTIqFwPBzHr1wPLq1cE0FwKhtPTg7qdPPx1tLaJAkAzLxUCA3nvW05aW4Nc4BJezzAUaR8gdCgTJGHc/MIaQa4GQrOeVV16JuJLsWrMmODq4shIGDoy6msBHPxpMn3wyCCyJjgJBMmbjxo289957lJeXM3z48KjLaeO0004D4O9//3vElWRX8jQR06dHW0eqKVOC7qsdO9RtFDUFgmRMrnYXAcyYMYOSkhJWrlxJQ0ND1OVkzauvBtNcCgQzuPzy4P5vfxttLb2dAkEyJpcDIRaLMXPmTNydl3rROZiTLYSTToq2jkN9/vPB9A9/gPffj7aW3kyBIBmTy4EAcOqppwK9p9vIPTe7jCA4avr004NTWPzhD1FX03spECQj3J0XX3wRgOm59u2T0NvGETZvhvp6GDwYKnLwMFF1G0VPgSAZUVNTw4YNGxg8eHDOBkKyhfDiiy/S2toacTWZlxw/OOmkoN8+13zmM1BcHAwsb90adTW9kwJBMmLRokUAnH/++TlzyopDjRo1inHjxrF7926qq3v+6bVytbsoacgQ+NjHIB6H+++PupreSYEgGZEMhAsuuCDiSg6vN3Ub5XogAFx/fTC9807YvTvaWnojBYKEbv/+/TzzzDMAfDR51FGO6o2BkGt7GKU6/3w47bRgrOOuu6KupvdRIEjonn/+efbs2cOUKVMYPbpbV1PNmt4SCA0NsH59cO6gXDiHUUfM4Kabgvt33AHvvRdtPb2NAkFCly/dRRDsEtuvXz/WrVvHunXroi4nY1asCKYnnghFoZz0PnPOOSc44d2uXfAf/xF1Nb2LAkFClwyEXO8uAigqKuKTn/wkAPfee2/E1WTOCy8E0xkzoq2jM8zgBz8I7t95J9TWRltPb6JAkFBt27aN5cuXE4vFOPPMM6Mup1O+8pWvAHDfffexf//+iKvJjOSZRM89N9o6Ouvss+Hii4Ourrlzobk56op6BwWChOrRRx8F4KyzzqKsrCziajrntNNOY/LkyWzfvp3HH3886nJC19gIzz0X/PLOl0AA+NWvYNQoeP55+O53o66md1AgSGj27NnDzTffDMBVV10VcTWdZ2YHWgnz58+PuJrwPfdccCH7GTOCff3zRXk5PPwwFBbC7bfDY49FXVHPp0CQ0PzkJz+hrq6Ok08+mTlz5kRdTpd84QtfIBaL8fTTT1NTUxN1OaFKdhclL0STT848ExK/MZgzB3pgXucUBYKEYseOHdx+++0A/OhHP8q7C9cPGjSIz33ucwDceeedEVcTrqeeCqZ5MMbfrn/+Z/iXf4HWVvjKV+Ab3wgu8iPhM3ePuoZ2VVVV+dKlSyNbvzu88w7U1AT7bzc2Bh/IoiIYOza44tTRR+f+LnzZcuONN3LXXXdx0UUXsXDhwqjL6ZZly5Zxyimn0Nrayp133snXv/71qEtKW10djBwJffoEn+fS0qgr6r5f/hKuvTa4qtrYsfD97wcnxNP/wbbMbJm7V3Xrze6e9g24EFgD1ADz2nm9FHgk8fpLwLgjLXPmzJmeTS0t7suXu999t/vcue4VFe5BLHR869fP/fzz3X/4Q/cXXwyW0dvs2rXLr7jiCgfczHzFihVRl5SWBx54wAEH/N57783KOuNx92XL3P/1X90//nH3KVPcBw50nzrV/frr3X//++5/tn7zm+CzeuGFoZYcmWefdT/hhIP/B0ePdr/uOvcnn3R/552oq8sNwFLv5nd52i0EMysE1gLnA7XAEuBSd1+VMs/1wFR3v9bM5gKfdPfPHW6506ad5N/85t2UlrYSi7VSWtpKuudIS/1b9+wpYt26/qxceRQrVx7FqlUDaWwsbjN/WVkLo0Y1MnJkE337tlBQ4DQ3F1BXV0ZdXRk7d8bazD9gwH5mzqzn5JN3MnPmTo46quv7yh367xGPw/btMTZt6kdtbV8aGop5//1i9u4NNoYZFBfHKStroayslbJCRBsLAAALWElEQVSyFvr0aSUWC6ZlZS3EYi2UlbVQWtpKa2sBLS0FNDcbzc0FNDcHj1Ofa2kpwB0KCpzCQk+ZQmGhA628884Otm6tZfHip9i+fTslJTG+9rUb+fjHZ7eJzXg8mCZrLSgIbsn7ZsE8yfkOvW8WDCoG6z54Px4PWmzJ+VNvEJw1s6jo4LSgINh1saMbQCwW/IJeuPC/uPvuO4B9HH/8OM4553SOPXYi+/aVsW9fKfv2xdi3r5SmpmIaG4toaiqkoMCJxYJ/h5Ejmxg1qomRI5uIxeLt/jvv32+sWDGIF14o54UXhrJjR6zd+ZImT36Xb31rNWPGNHbp83TbbZN5+umRXHvtWubM2dSl9+aq1lZ45pkR/OY3lWze3KfNa0OG7GXMmEaGDNnH4MH76dMn+O4oLY0fmBYVxTEDs+CDGdz/4GPwlNcOPk7OkzoNfkN88Pnk+9p/re3/9faW2dX1FRcb//iPp3S7hRBGIJwKfN/dL0g8/t8A7n5ryjyLEvO8YGZFQB1Q7odZ+dFHT/ZNm1Yd8mwT0AjsSbkd7vFeoAgoAQYAw4ARwPHAqHbWuh74W8rtDZIbvn0jgDOBswkaSeNTXosDrwEvEjSKVhM0kOo7WNZgYAwwAZgETExMjwf6dPAeyX2bgXVAA8HnsZTg3/gYgs9l0hbgCWAxwedkM8G//ZnADQSf1ybgW8A9nVx3IbAJGAmcCKxM70/JSTOB2cAFwBSgb7TlRK4BOKrbgRBG79togk9dUi1wSkfzuHuLmTUAQ4CdqTOZ2TXANQBDhowlFttKa2uMeDxGa2spUJa4pb/vXEHBXvr23chRR1Vz1FErGTiwmlgs9ct6dOJ2eGbvAo/h/hiNjRXs3Hky9fUns2vXVOLx6cB04NqU9e6nsLCRoqJG3AuIx0tpaelDPN5x525paT19+myib99NxGL1FBXtobBwL0FYGfF4CS0tMVpb+9DSUkZraxktLX1obY0lpmUH1mHWQkFBMwUFwfTg4+aUxy2JZRfgnrwVtpkWF5cRi/UhFuvL4MGDKSiwNr+mUlsAyVtqCyC1FVBY2HHrITlPa+vBW/I9ydbCoTf3oJ+5uTmYtrQE7ysu/uAt2YqAYKBy795gum8fNDXFeeedRhoa9hKP76e4eA+Fhcnb+xQU7Dnwb+FuiX/LfjQ1jaSpaTRNTSNw7/hz1K/fOoYO/TtDh/6d/v3fTPnFODhxA3iJ5uZq3nzzOurqLgTuprKyisrK3xzxs7l163msXj2SsrJaPvShYZidc8T35Ke/AH/B3WhqGkFT02j27x/E/v2DaW0tS3yHlCSmpcTjRUDyp7bhbu08bv81oM3rbR3u9a6+N3We9l5r//WSkv3s2kW3hREI7e1OcujP6s7Mg7vPB+ZDclB5ZMpr0NQUDO7u2XPwdujj1NvevcF/9pIS6NsXhg8PbsceC2PHxigoOA44Dvhkt//4jjQ1wbJlwSkDli4NBqdrauC990qIx0tobj6qzfwDBwZXsRo/HiZNCm4TJwa3o44aQhCCOXze4h6pAOiXuHVdSwts2gQbNgSfx6amIMQmTAg+g336HEPQUvhCp5b361/Dl74EGzZcyeWXX8kPftDxhW5aW2Hy5OD+T39awZVXLu7W3yD5J50d/MIIhFqCvo6kCoL2b3vz1Ca6jAYC73RlJWbBnhJ9+sDQoemUmx1lZXDGGcEtyT0IqffeC871XlQU9Fv36xfcpGcpKgr2RqusDGd5V1wR/MC5/HL44Q+DL/2bb27/C2DBAli7Nlj3ZZeFs37p+cIIhCXABDOrJOj4nAt8/pB5ngC+CLwAzAH+fLjxg57KLAiKsrKgpSLSVZdeGrQyPv95uOWWoBVy221tQyEeDwID4NvfPtglJnIkaR+Y5u4tBKNeiwhGThe4e7WZ3WRmsxOz/RIYYmY1wNeBeemuV6S3+uxn4ZFHghbI7bfDP/1T0HUKQRj8+7/D6tXBvvpXXBFtrZJfdGCaSJ56/PHgwvTNzcF5f669FhYtgpdfDl7/xS/g6qujrVGyL50D03TqCpE8dcklwXmKqqpgx46gm+jll4Mjk3/1K8ij8wtKjtBB3yJ57KyzghBYvBh++lOYOhW++U3tpCDdo0AQyXNmcN55wU0kHeoyEhERQIEgIiIJCgQREQEUCCIikqBAEBERQIEgIiIJCgQREQEUCCIikqBAEBERQIEgIiIJCgQREQEUCCIikqBAEBERQIEgIiIJCgQREQEUCCIikqBAEBERQIEgIiIJCgQREQHSDAQzG2xmfzKzNxPTQR3M90cze9fM/jud9YmISOak20KYByx29wnA4sTj9vwY+EKa6xIRkQxKNxAuAR5I3H8A+ER7M7n7YmB3musSEZEMKkrz/cPdfSuAu281s2HpLMzMrgGuSTzcZ2Yr06yvpxgK7Iy6iByhbXGQtsVB2hYHHd/dNx4xEMzsaWBEOy99p7sr7Yi7zwfmJ9a71N2rwl5HPtK2OEjb4iBti4O0LQ4ys6Xdfe8RA8HdzzvMireZ2chE62AksL27hYiISLTSHUN4Avhi4v4XgcfTXJ6IiEQk3UC4DTjfzN4Ezk88xsyqzOze5Exm9hzwe+BcM6s1sws6sez5adbWk2hbHKRtcZC2xUHaFgd1e1uYu4dZiIiI5CkdqSwiIoACQUREEiIPBDO70MzWmFmNmX3gSGczKzWzRxKvv2Rm47JfZXZ0Ylt83cxWmdlrZrbYzMZGUWc2HGlbpMw3x8zczHrsLoed2RZm9tnEZ6PazB7Kdo3Z0on/I0eb2TNm9mri/8nHoqgz08zsPjPb3tGxWhb4v4nt9JqZzejUgt09shtQCKwDxgMlwApg8iHzXA/8LHF/LvBIlDVHvC0+AvRJ3L+uN2+LxHz9gb8CLwJVUdcd4ediAvAqMCjxeFjUdUe4LeYD1yXuTwbeirruDG2LDwMzgJUdvP4x4EnAgA8BL3VmuVG3EGYBNe6+3t33Aw8TnA4jVerpMR4l2FPJslhjthxxW7j7M+7emHj4IlCR5RqzpTOfC4AfArcDe7NZXJZ1Zlt8GbjH3XcBuHtPPR6oM9vCgQGJ+wOBLVmsL2vc/a/AO4eZ5RLg1x54ETgqcazYYUUdCKOBTSmPaxPPtTuPu7cADcCQrFSXXZ3ZFqmuIvgF0BMdcVuY2UnAGHfv6WfQ7czn4jjgODN73sxeNLMLs1ZddnVmW3wfuNzMaoGFwNeyU1rO6er3CZD+uYzS1d4v/UP3g+3MPD1Bp/9OM7scqALOymhF0TnstjCzAuAnwJXZKihCnflcFBF0G51N0Gp8zsymuPu7Ga4t2zqzLS4F7nf3O83sVOA3iW0Rz3x5OaVb35tRtxBqgTEpjyv4YBPvwDxmVkTQDDxcUylfdWZbYGbnEZxHara778tSbdl2pG3RH5gCPGtmbxH0kT7RQweWO/t/5HF3b3b3DcAagoDoaTqzLa4CFgC4+wtAjODEd71Np75PDhV1ICwBJphZpZmVEAwaP3HIPKmnx5gD/NkToyY9zBG3RaKb5OcEYdBT+4nhCNvC3Rvcfai7j3P3cQTjKbPdvdsn9cphnfk/8hjBDgeY2VCCLqT1Wa0yOzqzLTYC5wKY2SSCQNiR1SpzwxPAFYm9jT4ENHjizNSHE2mXkbu3mNkNwCKCPQjuc/dqM7sJWOruTwC/JGj21RC0DOZGV3HmdHJb/BjoB/w+Ma6+0d1nR1Z0hnRyW/QKndwWi4CPmtkqoBX4lrvXR1d1ZnRyW3wD+IWZ/S+CLpIre+IPSDP7HUEX4dDEeMm/AsUA7v4zgvGTjwE1QCPwpU4ttwduKxER6Yaou4xERCRHKBBERARQIIiISIICQUREAAWCiIgkKBBERARQIIiISML/B6TNLZ+yHn73AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "t, C, X, u, LA = sp.symbols('t, C, X, u, LA')\n",
    "c = 0.3\n",
    "\n",
    "def init(x):\n",
    "    middle, width, height = 0.4, 0.1, 0.5   \n",
    "    return height/width**10 * (x%1-middle-width)**5 * \\\n",
    "                              (middle-x%1-width)**5 * (abs(x%1-middle)<=width)\n",
    "    \n",
    "dico = {\n",
    "    'box': {'x': [0., 1.], 'label': -1},\n",
    "    'space_step': 1./128,\n",
    "    'scheme_velocity': LA,\n",
    "    'schemes': [\n",
    "        {\n",
    "            'velocities': [1, 2],\n",
    "            'conserved_moments': u,\n",
    "            'polynomials': [1, LA*X],\n",
    "            'relaxation_parameters': [0., 2.],\n",
    "            'equilibrium': [u, C*u],\n",
    "            'source_terms': {u: -sp.Abs(X-t)**2*u},\n",
    "        },\n",
    "    ],\n",
    "    'init': {u: init},\n",
    "    'generator': 'cython',\n",
    "    'parameters': {LA: 1., C: c},\n",
    "}\n",
    "\n",
    "sol = pylbm.Simulation(dico) # build the simulation\n",
    "viewer = pylbm.viewer.matplotlib_viewer\n",
    "fig = viewer.Fig()\n",
    "ax = fig[0]\n",
    "ax.axis(0., 1., -.1, .6)\n",
    "x = sol.domain.x\n",
    "ax.plot(x, sol.m[u], width=2, color='k', label='initial')\n",
    "while sol.t < 1:\n",
    "    sol.one_time_step()\n",
    "ax.plot(x, sol.m[u], width=2, color='b', label=r'$D_1Q_2$')\n",
    "ax.legend()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
